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To advance hierarchial equations of motion as a standard theory for quantum dissipative dynam- 
ics, we put forward a mixed Heisenberg-Schrodinger scheme with block-matrix implementation on 
efficient evaluation of nonlinear optical response function. The new approach is also integrated with 
optimized hierarchical theory and numerical filtering algorithm. Difi'erent configurations of coherent 
two-dimensional spectroscopy of model excitonic dimer systems are investigated, with focus on the 
effects of intermolecular transfer coupling and bi-exciton interaction. 



I. INTRODUCTION 



Coherent multi-dimensional spectroscopies provide 
powerful tools to investigate various effects of molecu- 
lar interactions and dynamic correlations [ll-Q. Recent 
experiments show the evidence of long-lived quantum co- 
herence in photosynthesis systems 1-%, even at room 
temperatures [l^, EI | ■ Theoretical studies are also car- 
ried out towards the simulating and understanding of 
the involved excitation energy transfer processes ]\% - 
\vf\ . In such systems, the pigment-protein interaction 
has the same magnitude as the pigment-pigment transfer 
coupling, and the time scale of environment memory is 
about that of excitation energy transfer. These charac- 
teristics require the non-Markovian and non-perturbative 
quantum dissipation methods. Among them^ hierarchi- 
cal equations of motion (HEOM) approach |18l - [23| , an 
equivalence to the Feymann- Vernon influence functional 
path integral (23 - [26| but numerically more efficient al- 
ternative, emerges as a standard method. 



There are many recent efforts, in both theoretical and 
numerical aspects, on advancing exact HEOM to be a 
power and versatile tool for the study of various quan- 
tum dissipative systems. HEOM formalism was origi- 
nally proposed in 1989 by Tanimura and Kubo for semi- 
classical dissipation 27^ . Formally exact HEOM formal- 
ism (l8l - l23i] for Gaussian dissipation in general, including 
its second quantization [2^, has now been well estab- 
lished. The major obstacle of HEOM is its numerical 
tractability. The number of equations involved in the 
theory are usually huge. Brute-force implementation is 
greatly limited by both memory and processing capabil- 
ity. One major numerical advancement is the on-the-fly 
filtering algorithm [29}. It goes with a preselected error 
tolerance on the properly scaled HEOM. The filtering 
algorithm dramatically reduces the effective number of 
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equations, and it by nature also automatically truncates 
the hierarchy level. In the formulation front, the Fade 
spectrum decomposition (FSD) scheme j30, 31] has been 
proposed for an optimized HEOM construction [S^]. It 
dramatically reduces the number of equations, in com- 
parison with that of the conventional Matsubara expan- 
sion based formalism at same accuracy. A priori accuracy 
control criterion has also been proposed |32h34I 1 so that 
the optimized HEOM can be used confidently, without 
costly convergency test. 

In this work, we discuss two additional techniques 
to further improve the efficiency of HEOM in evaluat- 
ing such as the third-order optical response functions 
and two-dimensional spectroscopy. One is the mixed 
Heisenberg-Schrodinger scheme. The conventional ap- 
proach follows the direct realization of third-order re- 
sponse function from the view of Schrodinger picture. It 
propagates the reduced system density operator in three 
nested time intervals, denoted as ti, t2, and is, between 
three interrogations and the time of detection which is 
through dynamic variable such as transition dipole. In 
the new scheme, while retaining ti- and i2-pi'opagations 
on state variables, we unlash from the nested loop via 
its propagation in Heisenberg picture on dynamic vari- 
ables. HEOM in Schrodinger picture and Heisenberg pic- 
ture are detailed in Sec. Ill Al and Sec. Ill B| respectively. 

Another advancement in this work is the block-HEOM 
dynamics, for its efficient evaluation of coherent two- 
dimensional spectroscopy, see Sec. lIIII In optical pro- 
cesses, the system Hamiltonian is considered to be block 
diagonalized in electronic manifolds, in virtue of Born- 
Oppenheimer principle. We also assume that the re- 
laxation between different manifolds is negligible. The 
resulting HEOM dynamics will be in the block-matrix 
form, even as they involve in the third-order optical re- 
sponse functions, where the dynamics occur in popula- 
tion and/or coherence states in electronic space. Combin- 
ing the block-HEOM theory with the mixed Heisenberg- 
Schrodinger scheme greatly facilitates the evaluation of 
third-order optical response functions. Moreover, numer- 
ical filtering algorithm [29j is now also extended to the 



2 



present formalism. 

We exemplify our method with model excitonic dimer 
systems in Sec. lIVI The well-established ki, kn and kni 
types of coherent two-dimensional spectroscopy [5] are 
evaluated with the present HEOM dynamics that inte- 
grates all the state-of-the-art techniques. Finally we sum- 
marize the paper in Sec.lVl 



II. HIERARCHICAL EQUATIONS OF MOTION 

A. HEOM in Schrodinger picture 

HEOM couples the reduced system density operator of 
primary interest, p{t) = trBPtotai(i), to a set of auxiliary 
density operators (ADOs). As an exact and nonpertur- 
bative theory, HEOM accounts for the combined effects 
of system-bath coupling strength, environment memory 
timescales, and many-body interactions [l^, [2^. The ex- 
plicit form of HEOM is defined upon a certain statistical 
environment bath "basis set" that decomposes the inter- 
acting bath correlation functions into distinct memory- 
frequency components. Without loss of generality, we 
exemplify it with the case of single-mode system-bath in- 
teraction, H'{t) — —QF^{t), where Q and F^{t) are oper- 
ators in the reduced system and the stochastic bath sub- 
spaces, respectively. In general, H'{t) can be expressed 
in multiple-modes decomposition form. The stochastic 
bath operator F^{t) assumes a Gaussian process. The 
influence of bath is therefore described via its correla- 
tion function C{t) = {Fji{t)F^{Q))^. It is related to the 
bath spectral density J{uj) via the fluctuation-dissipation 
theorem [25,. ,35J: 



1 

C{t) = - / 
J-, 



did- 



(1) 



where 1/(1 — e~"") is Bose function at inverse tempera- 
ture P — h/{kBT). We set ?i = 1 hereafter. 

To construct a HEOM formalism, one need to expand 
C{t) in a finite exponential series, on the basis of certain 
sum-over-poles scheme, together with the Cauchy residue 
theorem of contour integration applied to Eq. ([T]). In this 
work we adopt the Drude model. 



J{u) 



(2) 



In this case, [N/N] Fade spectrum decomposition (PSD) 
for Bose function is shown to be the best (30l - l34| . It 
results in an exponential expansion of interacting bath 
correlation function 1321 ■ 



N 



C(0 «^Cfce-^'=*+2AA„5(t), 



fc=0 



with 



A 



N 



A^7 



2{N+l){2N + iy 



(3) 



(4) 



The fc = term with 70 = 7 is the Drude pole contri- 
bution, while other N contributions with k = 1, • ■ • ,iV 
are from the [N/N] PSD Bose function poles jjifc/o} that 
are all positive and can be easily identified [31'|. The 5- 
function term in Eq. ([3]) re-sums the off-basis-set residue 
outside the finite sum-over-poles scheme. This is the only 
approximation involved not just in the bath correlation 
function, but also in the resulting HEOM that goes there- 
fore by a convenient accuracy control criterion prior to 
dynamics evaluation [b^ - Is^ . 

The exponential expansion form of C{t) dictates the 
construction of HEOM. According to Eq. ([3|) it reads ex- 
plicitly as m,!!!,!!!!!!] 



p„=-[iC{t)^-1„+ ^NQ^]Pn 
N 

N 

-i^ VK + l)|c/c|QPn+' 

fc=0 



(5) 



with Cpn = [H, pn] for the reduced system Liouvillian, 
Qpn = [Q,Pn], and 



7n 



N 
k=0 



(6) 



The ADO labeling index is now specified as n = 
{riQ, • • • , uk}, with rifc > 0. It consists a set of non- 
negative indices, following the exponential expansion of 
Eq. The index in the last two terms of Eq. ([5]) dif- 
fers from n only by changing the specified Uk to ± 1. 
Let no + rti + • • ■ + uk — n and call the individual 
Pn = Pno,ni,--- .UK ^u n"^-tier ADO. It depends on its 
associated {n ± l)"^-tier ADOs, specified individually by 
the last two terms in Eq. ([S]). The reduced system density 
operator of primary interest is just the zeroth-tier ADO, 
po(i) = p{t). The AAT-term in Eq. (O arises from the 
white-noise residue in Eq. All ADOs in Eq. ([5]) are 
dimensionless and scaled properly to support the efficient 
HEOM propagator via the on-the-fly filtering algorithm 
that also automatically truncates the level of hierarchy 
[29| . The nonperturbative nature of HEOM has been dis- 
cussed in detail, see Ref. for example. HEOM ([5]) is 
called to be in Schrodinger picture, for its governing the 
system state variables, which are the ADOs including the 
reduced system density operator of primary interest. 



B. HEOM in Heisenberg picture 

HEOM supports the evaluations of not only the ex- 
pectation values but also the correlation and response 
functions of any dynamics variables of the reduced sys- 
tem. Nonlinear optical response functions for two- 
dimensional spectroscopy will be described in the next 
section. They will be evaluated based on the mixed 
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Heisenberg-Schrodinger scheme and the block-HEOM 
dynamics. 

Note that HEOM consists a set of hnear coupled equa- 
tions. The HEOM space algebra is mathematically the 
same of linear space, regardless its physics contents. We 
recast HEOM as p{t) = —A{t)p{t), with the column 
vector p = {pn=o,pn^o} of ADOs defining a basic ele- 
ment in HEOM space. The dynamic generator A{t) is 
determined by specific form of HEOM, i.e. Eq. The 
HEOM propagator Q{t,T) by which p{t) = g{t,T)p{T) 
satisfies dQ{t,T)/dt = —A{t)Q{t,T). In the absence of 
time-dependent external field, jC,{t) = £ and therefore 
A(t) = A are time-independent, resulting in Q{t,T) = 
exp[-A{t - t)] = g{t - r). 

The expectation value of a system dynamical variable 
reads 

A = tr(^p) ^ = {{A\p)) ^ ^((^nlPn)). (7) 

all n 

The first two identities are the conventional reduced Li- 
ouville space expressions. The last two identities are the 
extensions to HEOM space, in which, according to the 
third identity. 



A = {A=o = A, An^o = 0} = A{0). 



(8) 



The last identity will be used later to specify the initial 
values of HEOM in Heisenberg picture. 

HEOM supports also the evaluation of correlation and 
response functions. For example, we apply linear re- 
sponse theory to HEOM, resulting in linear correlation 
function the expression of 



CAB{t) = {{A\g{t)B\p°^ 



(9) 



which is equivalent to CAsit) = {{A\QM{t)\Bp''j^j/ 
{A{t)B{0)). The latter is the conventional expression, 
defined in the full system-plus-bath material space. How- 
ever, Eq. (jni) is defined in HEOM space of the reduced 



system only. The equilibrium p°'^ = {Pp^O'/^n%} 
the stationary solution to HEOM, i.e. Apoq = 0, to- 
gether with the normalization condition of the primary 
reduced system density operator Trpo = 1- The result- 
ing {Pp^o 7^ 0} in general account for the initial system- 
bath correlations. In Eq. dH), = {BPn=0' 
in comparison with the full material space counterpart 

of BpXf = -^Pm- Denote also Bp = {pn=oB, Pn^oB} for 
later use; cf. Eq. p^ . 

Apparently, the propagator Q(t) can take action on a 

system state variable, e.g., Q{t){Bp'^^) in Eq. ©. This 
is the Schrodinger picture. It can also act from right 
to left on a dynamic variable, i.e., A(t) = AQ{t). This 
is the Heisenberg picture, the HEOM analogue of the 
conventional A{t) = AQM{t) that satisfies the Heisenberg 
equation, A = —iACm = — iJ^f], defined in the full 
system-plus-bath material space. 



HEOM in Heisenberg picture satisfies A{t) = —A{t)A, 
with A{0) — A defined in Eq. ([5]). Its explicit expressions 
can be obtained as follows. Let us start with the identity. 



:(A|A|p)) = {{A\p)) = {{A\p)), 



(10) 



where A = AA and p = Ap, with A acting from right 
and left, respectively. Note that A = A and p = p. 
From Eq. ([5]) we have 



all n 
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fc=0 
N 



lJ2V{nk + l)\ck\{{A,\Q\p,p 



k=0 



(11) 



In contact with the second quantity in Eq. (|10p , we recast 
every individual term above with respect to the same 
Pn = p„o,ni- - by using for example the identity. 



We can therefore recast Eq. (|lip as 



A\p))=J2{U'^\'^ + l'^+^NQ^\Pn) 
all n 



N 



fc=0 ' 
N 



i rik + l 
|cfe| 



(cfc(M„+Q|pn))-4((Q^,+ |pn)) 



+ iJ2^''k\ck\ {{A-\Q\p„) 



k=0 



Accordingly, HEOM in Heisenberg picture reads explic- 
itly as [37|: 



fc=0 
N 



|Cfe| 



-iJ^^y^klcklA-Q. 



(12) 



k=0 



Here, OC = [O, H] and OQ = [O, Q], following the iden- 
tities of CO = [H, O] and QO = [Q, O] defined earlier. 
Apparently, Eq. ([SJ and Eq. (IT^ are equivalent but just 
in different pictures. 
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FIG. 1: Eight Liouville-space pathways (upper) and double- 
sided Feynman diagrams (lower) for two-dimensional spec- 
troscopy in the rotating wave approximation. Each Liouville- 
space pathway starts from the upper-left circle, while the 
double-sided Feynman diagram starts from the bottom, fol- 
lowing the convention of Ref. yj. 



III. EFFICIENT HEOM EVALUATION OF 
TWO-DIMENSIONAL SPECTROSCOPIES 

A. Nonlinear optical response functions via 
block-HEOM dynamics 

Now turn to the third-order optical response functions, 
as probed by coherent two-dimensional spectroscopies, 
operated with short pulsed fields in certain four-wave- 
mixing configurations 0, Q . Following the similar alge- 
bra of the linear correlation/response function, the third- 
order optical response function in HEOM space is ob- 
tained to be 

R^'Ht,,t2,h) = {{iJi^^\g{h)v^,g{t2)v^,G{t{)'D^,\p-^)). 

It is just the HEOM space analogue of the conventional 
full system-plus- bath material space expression Ap- 
parently the mixed Heisenberg-Schrodinger scheme as 
commented in Sec. U will greatly facilitate the evaluation 
of third-order response function. 

As a kind of four-wave-mixing spectroscopy, the 
wavevector ks of the signal field satisfies the phase- 



matching condition. It is that = itka ± k2 ± ki, 
in relation to the three incident pulsed fields interact- 
ing with the system sequentially. Three basic configura- 
tions of coherent two-dimensional spectroscopy 0, 1^ will 
be classified in Sec. lIIIB] Their efficient evaluation via 
block-HEOM in mixed Heisenberg-Schrodinger scheme 
will be detailed in Sec. lIII Cl 

In the following, we show that the optical response 
function can be recast in the block-HEOM dynamics 
form. Consider the third-order optical processes involv- 
ing the initial ground \g), the excited |e), and doubly- 
excited I/) manifolds of electronic states. Assume also 
that the relaxation between different manifolds is negli- 
gible. This implies that not only the system Hamiltonian 
but also the dissipative mode Q are block diagonalized, 
in virtue of Born-Oppenheimer principle. 

On the other hand, the transition dipole operators in- 
volved in the third-order optical response function are 
also in the block-matrix form, although not diagonal. For 
example, the transition dipole /x_ = jlge\g){e\ +/ie/|e)(/| 
amounts explicitly to 



^ffe Og/ 

Oeg Oee Ae/ 
.6/9 6/e 6//. 



The matrix Ouv h^ the (uf)-block, with u,v G {5,6,/}, 
has the order of NuX Ny. We have Ng = 1, Nf, = M and 
Nf — M{M — l)/2, for a molecular aggregate of size M, 
assuming the simple exciton model. 

Expanding the third-order optical response function 
in the eight Liouville-space pathways, as shown in FigJT[ 
and their complex conjugate counterparts leads to [ll l38l| 



R^^Kh, t2, h) = i^Y. iRc.{h,t2,h) - c.c.]. (13) 

These eight pathways contributions are expressed in 
terms of the block-matrix dynamics in HEOM space as 



Riit-s, t2 


ti) = 


{{f^ge\Geg{h)^legGee{t2)^lgeGeg{tl)'jl,g\pll))e- 




^2(^3, h 


ti) = 


{{t^ge\Qeg{t3)'fl^ggee{t2)~i^egGge{h)p,gJpll))e- 




Rsih, h 




{{fJ'gelQeg (^3 ) /^eg ^33 (^2 ) Meg ^se 1 ) Mge 1 Pgg )) ^ 


-icJcg(t3— tl) 

f 


R^its, h 




{{tJ'ge\0eg{t3}]^egSgg{h)]^ge0eg{tl)l^eg\Pg1))^' 


-icJeg(i3+*l) 


Rbits, h 


ti) = 


-{{^■ef\0fe{h)p-f^Qee{t2)'P'egQge{tl)V'ge\p7g)) 




Reit^, ^2 


ti) = 


-{{f^ef\Ofe{h)'flf,gee{t2)^lge0eg{tl)'fleg\p7a)) 




Rrits, h 


ti) = 


-{{Pef\Q le{h)PgeQ Sa{t2)'tif,Qeg{tl)^eg\p7a)) 




Rsita, t2 




{{P'ge\Qeg{H)tefQf9ii2)PfeQeg{tl)ll^g\p7g))e 





(14) 
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The underlying optical processes will be discussed in 
Sec. lIII B1 The electronic phase factor in each indi- 
vidual Ra can be formally absorbed into the involving 
Green's functions, i.e., ^„„(i)e~*"""* — >■ with 
ujuv = — denoting the chosen reference frequency 
for the optical transition between two specified electronic 
manifolds. In the present notion, Quv{i) involves only 
slow motion dynamics, as the highly oscillatory optical 
frequency component is factorized out for numerical ad- 
vantage. The corresponding block HEOM dynamics in 
both Schrodinger and Heisenberg pictures will be detailed 
in Sec. lIII Cl There are other advantages for the present 
notion. The overall electronic phase factor can in fact be 
used to distinguish rephasing versus non-rephasing opti- 
cal processes [ll, [S^, [S^ and also to visualize the rotating 
wave approximation as seen below. 



B. Coherent two-dimensional spectroscopies 

There are three basic configurations of coherent two- 
dimensional spectroscopy, and their signals are denoted 
as iSki , 5'kii , and , respectively Q • For simplicity we 
adopt the rotating-wave approximation and the impul- 
sive fields limit. 

The ^ki signal goes with = ka + k2 — ki , the stimu- 
lated photon echo or rephasing configuration 39], while 
the S'kjj signal goes with k^, = ka — k2 -I- ki and is non- 
rephasing. With the aid of the double-sided Feynman 
diagrams in Fig.[l] these two signals are identified to be 

/>oo />oo 

5k,/„(c.3,i2,^i) = Rc / dtJ dtie'^^^''-'^'^'''^ 
Jo Jo 

xi?k,/„(i3,i2,ii), (15) 

and related respectively to 

Rki = R2 + R3 + R5 (rephasing) , ^^^^ 

^kii = -Ri + -R4 + i?6 (non-rephasing) . 

The rephasing versus non-rephasing nature of individual 
Ra in Eq. or Eq. ([T5| can be inferred easily from 
its overall electronic phase factor [l|, HI, . The signs 
associating with the frequencies 0^3 and uji in Eq. (jl5p 
are resulted from the incident fields in the specified four 
wave mixing configuration in the impulsive limit, as im- 
plied in Fig.[TJ These signs are just opposite to those 
in the electronic phase factor of participating Ra contri- 
butions. Thus, the participated pathway contributions 
via the electronic rotating wave approximation are also 
evident in Eq. (fTS)) . 

Experiments can also be performed in the configura- 
tion that the pulsed k2 -field is applied continuously not 
only after but also before the ki-field. The resulting 
signal amounts to 5ki+kii ~ S'ki + S'kn. It is in fact 
the pump-probe absorption configuration, involving all 
the six pathways Ri to Rq contributions. As inferred 
from Eq. (ITil) , these six pathways can be classified into 



the excited-state emission (i?i,i?2), ground-state bleach- 
ing (i?3,i?4), and excited-state absorption (i?5,i?6) con- 
tributions. In fact, the ti and t^ represent the excitation 
and detection time periods, and therefore, the wi and ^3 
in Eq. (jl5l) are the excitation and detection frequencies, 
respectively. The <2 denotes the waiting time, during 
which the system is either in the excited or the ground 
state manifold, with underlying dynamics being governed 
by Qee{t2) or 0gg{t2), respectively; see Eq. (jl4[) . 

The S'kiii signal goes with ks = — k3 -I- k2 + ki, the 
double-excitation configuration, and is related to 

Rkui = R7 + Rs, (17) 

as inferred from the double-sided Feynman diagrams in 
Fig.[T] The i?7 and Rs are the double- excitation absorp- 
tion pathways, involving the |/) |e) <— \g) processes, 
while the bra state remains in {g\. During the time ^2 
period, the km configuration explores therefore double 
quantum coherence dynamics governed by Qfg{t2), in 
contrast to the kj/n scheme involving electronic state 
population dynamics. The detection k3 field involves 
single-excitation absorption of (e| ^ {g\ in Rj and single- 
excitation emission |e) ^ \f) in _Rg, as evident in the 
corresponding double-sided Feynman diagrams in Fig.[T] 
The above analysis justifies the fact that the km-signal 
is designed to probe the correlation between single and 
double excitations [I^llSim. The two-dimensional half- 
Fourier transforms are therefore performed with t2 to re- 
solve the double-excitation frequency and with either t^ 
or ti to resolve the specified single-excitation frequency. 
In this work, we choose 

/•oo poo 

5k,„(c^3,c^2,ii)=Re / dtJ di2e^("^*^+'"^'^) 
Jo Jo 

X i?k„:(<3,i2,il). (18) 

In the independent exciton limit, the kni-signal vanishes. 
This can be seen from the involving i?7 and Rg contribu- 
tions [cf. Eq. (ITi)) ]. Besides the signs, these two contribu- 
tions differ by their single coherence dynamics in the t^ 
period, which is Qfe{tz) in i?7 but Q^g(t^) in Rg. These 
two contributions would cancel each other in the ab- 
sence of both inter-exciton transfer coupling and double- 
exciton correlation. Thus the kni technique serves as a 
sensitive probe for interactions between excitons. 

C. Implementation with block- HEOM in mixed 
Heisenberg— Schrodinger picture 

To implement the third order optical response func- 
tions in Eq. (jl4p . we start with the thermal equilibrium 
Pg^ in the ground-state |(7)-manifold. As described ear- 
lier, it is determined by the steady state solution to 
HEOM, involving now only the {gg)-hlock part. For the 
simple exciton model, the |(7)-manifold contains only one 
level, and the ADOs in {gg)-h\ock are all 1 x 1 matrices, 
resulting in = {pS"* = 1' Pn^'r = 0}- 
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Block-matrix multiplications are then followed: 

Pe<,(0). 



We can recast Eq. ([14]) as (up to the phase factors) 



l^egPgg — \H-egHn=0 'h^egPn^O I 



(19) 



Denote also Pge(O) = MgePgg = {Pn=o''Ase,Pn^o'^Ase}- 

They are the initial states for the block-HEOM Quv{ti) 
propagations in Eq. (IT4l) . Each ADO in puv is an x Ny 
matrix, as inferred in the (uw)-indexes, and also Hermite 
conjugate with its counterpart in i.e., /5„„ = {Pn"} = 

The ti- and <2-pi'opagations in Eq. (|14l) are imple- 
mented in a nested manner in Schrodinger picture. This 
picture is defined via the action- from- left of Quv{t) on 
state variables; e.g., Puv{t) ~ Guv{t)puv{0)- The block- 
HEOM in Schrodinger picture can be reduced from 
Eq. ^ as 



p7 



N 




iy^^{nk + l)\ck\QuvPll 

1 ^ 1, 



(20) 



fc=0 



Here, Huv^uv — ^uu^uv O^v^iw and Quv^uv — 
Quu^uv ^uyQvv Equivalently, O^^C^y^^ — OyyH^y 
HyyOyu and OyuQuv = OyuQuu " QwOyy, wHch wlU bc 
used in the following Heisenberg picture. 

The t3-propagation in Eq. ()14p is implemented in 
Heisenberg picture, in parallel with the ti- and t2- 
propagations. The Heisenberg picture is defined via the 
action- from-right of Guv{t) on dynamic variables, i.e., 
Ai,u{t) — AyyQyyit). Thc iultial condition is Ayu({)) = 
Ayy - {A™o = = 0}, following Eq.®. In 

consistent with Eq. or Eq. dll]), the block-HEOM in 
Heisenberg picture reads 



Ar=-*Ar(/:™+7n+A^QL) 
— i 

k=Q 
N 

fc=0 



|Cfc| 



(21) 



To evaluate the third-order optical response func- 
tions in Eq. ([T4|) with block-HEOM in mixed Heisenberg- 
Schrodinger scheme, we introduce 



Puv{t2',ti) = Guv{t2)Puv{0',ti) 

for three types of initial t2 conditions: 

Pee(0;tl) = PgePegih), 
Pgg{0;ti) ^ Pg^pegih), 
PfgiO;tl) = Pf^Pegih). 



(22) 



(23) 





i2iii/ 


= {{Pge{h)\P'egPee{t2]ti))) , 




t2,ti, 


= {{P-ge{h)\'p.egPL{t2]ti))), 


Rsits 


t2,ti) 


= {{tJ^ge{h)\piegp\g{t2-M))), 


Ri{t^ 


t2, til 


= {{Pg^{t^)\p^gPgg{t2]tl))) , 


i?5(i3 


t2,tif 


= -{{Pcf{h)\P fePle{t2]ti))) , 


Reih 


t2,ti) 


= -{{Pef{h)\P fePee{t2]ti))) , 


Riih 


t2,ti) 


= -{{Pef{h)\'PgePfg{h]ti))), 


Rsih 


t2,tij 


= {{P-ge{h)\p.efPfg{t2]ti))). 



(24) 



These are the final expressions for the mixed Heisenberg- 
Schrodinger scheme block-HEOM evaluation of third- 
order optical response functions and coherent two- 
dimensional spectrums via such as Eqs. ([T5l) and ([T8)) . 

We have also successfully extended the on-the-fly fil- 
tering algorithm to the block-HEOM dynamics in 
Eg. There involve the nested (^2; ii)-propagation 

in Schrodinger picture [Eq. (|20)) ] and the separated ^3- 
propagation in Heisenberg picture [Eq. ([21])]. Setting 
the filtering error tolerance at 2 x 10~^ is found to be 
sufficient for HEOM dynamics in both Schrodinger and 
Heisenberg pictures, as tested extensively on various sys- 
tems, with numerical accuracy by eyes. 



IV. NUMERICAL DEMONSTRATIONS 

We exemplify the efficient evaluation of coherent two- 
dimensional spectrums with a model excitonic dimer. Its 
Hamiltonian reads 

H = eih\bi + €2blb2 + V{b\b2 + blh) + Ub\hblb2, 

where 6j„ (bm) denotes the exciton creation (annihila- 
tion) operator on the specified molecular site. The model 
system consists of a total four levels: \g) = |00) in 
the ground-state manifold, |e) — [10} and |01) in the 
single-exciton manifold, and [/) — |11) in the double- 
exciton manifold. The Rabi frequency within the single- 
exciton manifold is v^ei--e2p+4V^- The electronic 
transition dipoles fieg — (a^i|10) + /i2|01))(00| and pfe — 
|ll)((10|//2 + (01|/^i). Co-linear (xxxx) field polarization 
configuration is adopted, so that the effect of dipole di- 
rections on spectroscopic signals can be neglected. 

Each on-site transition energy experiences fluctua- 
tions, brought in through Qm = ^m^m influence of bath 
in Drude model. We neglect the cross correlation between 
different on-site fluctuations, and also the static disorders 
that are irrelevant to the methodology of this work. In 
the following, we set A = 60 cm^^ and 7^^ = 100 fs for 
each individual on-site Drude dissipation [Eq. ([2])] and 
77 K for temperature. In all cases, the [1/lJ-PSD scheme 
is sufficient according to the established accuracy control 
criterion fH-lll. 



7 




FIG. 2: Coherent two-dimensional spectra for the dimer sys- 
tem: ei — €2 — e, V = — 300cm~^, and U — 0, with 
Ml/pa ~ —5, at temperature 77 K. Drude dissipation pa- 
rameters for each on-site excitation energy fluctuation are 
A = 60cm~^ and 7"^ — 100 fs. Frequencies (when appli- 
cable) are reported in terms of Aoji — ui — e, Alj^ — lus — e 
and Alu2 = uj2 — 2e, as the electronic reference transition fre- 
quencies are chosen to be uieg = ojfe ~ ^ and LOjg = 2t. 



(a) S,(co3,t,,co, 




-400 400 

Aw (cm~^) 
(b) S,_(a>3,t2,(0,) 

t, = 










Ss"""^ 




-800-400 400 800 

A(o (cm" 



-800-400 400 800 

Am (cm"^ ) 




-1600-800 800 16D0 -1600-800 800 1600 -1600-800 800 1600 

Aw^(cm"^) A(o^(cm"^) Aw^(cm"^) 

FIG. 3: Same as Fig.[2]but with a finite bi-exciton interaction, 
U = 200cm"\ 



FIG. 4: Same as Fig.[3]but for ex = e - 300cm"\ £2 = e + 
300 cm~^, and V = Q, with ^i//i2 ~ 1. 



Figures [2H1] exemplify three representative cases of 
^ 0, C/ = 0), {V ^Q,U ^ 0), and {V = Q,U ^ 0), 
respectively, but sharing a common value of Rabi fre- 
quency (ei — £2)^ + ^V^. Thus the peaks and valleys, 
which reflect the transition frequencies between nonlocal 
eigenstates, are distributed at similar positions in these 
three figures. Highlighted are therefore the effects of ex- 
citonic transfer coupling V and bi-exciton interaction U 
on 5'ki, 'S'kii, and 5*1:1115 as depicted in the panels (a), (b) 
and (c) of each individual figure. Intensities are reported 
in their relative values, with a common factor over all 
frames in these figures. 

General speaking, both V and U affect correlations be- 
tween different monomers, manifested via the cross peaks 
in S'ki and S'kii- ^he other hand, Skm that vanishes 
when U = V = is specialized for the correlation be- 
tween single and double excitation coherence. It provides 
a different way to visualize the effects of finite U and/or 
as evident from Figs.IlHH The presence of U mani- 
fests mainly on the separation of peaks and valleys (cf. 
Fig.Hlfor example) and it is particularly prominent in the 
Skin spectrum. 

Other observations in rephasing S'ki and non-rephasing 
Skii are briefed as follows: (i) The peaks arise from ex- 
cited state emission {R2 and Ri) and ground state bleach- 
ing (i?3 and i?4), while the valleys are from excited state 
absorption (i?5 and Rq); (ii) The peaks/valleys are along 
the diagonal direction in Ski and anti-diagonal in Sk„ , in 
line with Eq. (jl5p . Inhomogeneity that is not included in 
calculations affects mainly the diagonal direction. Thus, 
it elongates the diagonal peaks in Ski ; while smears those 
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in 'S'kii; Bi-exciton interaction U shifts the excited 
state absorption, with respect to the excited state emis- 
sion and ground state bleaching. These two components 
differ in signs, being of valley versus peak. Thus, can- 
celation would occur at least partially when U — Q] (iv) 
Excitation energy transfer is observed as the evolution of 
peaks/valleys intensities. This process is mainly respon- 
sible hy V 0; {v) The correlation effects arising from U 
are separated out in Fig.|4l The negative peaks in Fig.|4] 
(a) and (b) would cancel completely with the positive 
ones when U ~ 0. The larger the U is, the bluer shift 
of the negative peaks from their positive counterparts. 
Note that the dimer system studied in Fig.|4] is nonde- 
generate, rather than the degenerate ones in Figs. [5] and 
121 for the appearance of correlation and coherence be- 
tween two distinct monomers. 



V. CONCLUDING REMARKS 

In summary, we propose a mixed Heisenberg- 
Schrodinger scheme and block-HEOM theory, and 



demonstrate it with efficient evaluation of third-order 
optical response function and coherent two-dimensional 
spectroscopy. The new development has also been in- 
tegrated with the efficient numerical filtering algorithm 
j29| and the optimized hierarchical theory ^32i-i34i] . This 
is the state-of-the-art HEOM approach. For example, 
the calculations of all frames in Fig. [5] take only about 
three minutes of CPU time on a single processor of In- 
tel(R) Xeon(R) E5472 (3GHz). The development made 
in work will greatly facilitate the use of HEOM, an ex- 
act and nonperbative quantum dissipation theory, to the 
study of realistic systems. 
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